        SUBROUTINE CalDdatSigma(dall,obst,cbst,sigmaT,meandeltaT,stddeltaT)
        IMPLICIT NONE
 ! INPUT
        INTEGER dall
        REAL obst(dall),cbst(dall)
! OUTPUT
        REAL  sigmaT(dall)
        REAL meandeltaT,stddeltaT
! PARAMETER
        INTEGER jj,kk,ii
        INTEGER i,j,k
        REAL twostdratio
        REAL :: deltaT(dall)
        deltaT=0
        meandeltaT=0

        DO i = 1,dall
            deltaT(i)=abs(cbst(i)/obst(i))
            meandeltaT=meandeltaT+deltaT(i)
	ENDDO
        meandeltaT=meandeltaT/dall
	stddeltaT=0
        DO i=1,dall
          stddeltaT=stddeltaT+(deltaT(i)-meandeltaT)**2
        ENDDO
	 stddeltaT=sqrt(stddeltaT/dall)
	 DO i=1,dall
	       twostdratio=abs(deltaT(i)/(1.5*stddeltaT))
	        IF (twostdratio.GT.1.0)THEN
	            sigmaT(i)=stddeltaT*obst(i)*exp(twostdratio-1)
	        ELSE
	             sigmaT(i)=stddeltaT*obst(i)
	         ENDIF
          ENDDO
          END SUBROUTINE


! use- by subroutine: FwdInvResData
        SUBROUTINE CalGcsReslNorm(maxvp,dall,GGc,GGs,dv,sigmaT,Tdata,fwdTaa,resbst,res2Norm,resNSigma2Norm,PreRes)
        use  lsmrblasInterface, only : dnrm2
        IMPLICIT NONE
! INPUT
        INTEGER,INTENT(IN):: dall,maxvp
        REAL,INTENT(IN):: GGc(dall,maxvp),GGs(dall,maxvp)
        REAL,INTENT(IN):: dv(maxvp*3),sigmaT(dall),Tdata(dall)
! OUTPUT
        REAL:: resbst(dall)
        REAL:: res2Norm,resNSigma2Norm
        REAL:: PreRes
! PARAMETER
        REAL VsInv(maxvp),GcInv(maxvp),GsInv(maxvp)
        REAL fwdTgc(dall),fwdTgs(dall), fwdTvs(dall),fwdt(dall)
        REAL fwdTaa(dall),res(dall),resW(dall),per_res(dall)
        INTEGER jj,kk,ii
        INTEGER i,j,k
!--------------------------------------------------------------------!
! Target:
!  Calculate traveltime resiudal between observe traveltime
!  and forward traveltime which is generated by inversion result
!  using G from ref  model.
!   r=Gm-d
!--------------------------------------------------------------------!
        GcInv=dv(1:maxvp)
        GsInv=dv(1+maxvp:2*maxvp)
!--------------------------------------------------------------------!
! forward residual
        fwdTgc=MATMUL(GGc,GcInv)
        fwdTgs=MATMUL(GGs,GsInv)
        fwdTaa=0
        res=0
        resW=0
        per_res=0
        resbst=0
        DO i=1,dall
                fwdTaa(i)=fwdTgs(i)+fwdTgc(i)
                resbst(i)=Tdata(i)-fwdTaa(i)
                resW(i)=resbst(i)*1/sigmaT(i) ! here, 1/sigma is for
                per_res(i)=resbst(i)/Tdata(i)
        ENDDO
        res2Norm=dnrm2(dall,resW,1)
        resNSigma2Norm=dnrm2(dall,resbst,1)
        PreRes=sum(abs(per_res(1:dall)))/dall
        write(6,*)' Mean Residual Percentage (%)',PreRes*100
        write(66,*)' Mean Residual Percentage subroutine (%)',PreRes*100

        END SUBROUTINE

        SUBROUTINE CalVsReslNorm(maxvp,dall,GVs,dv,sigmaT,Tdata,fwdTvs,resbst,res2Norm,resNSigma2Norm,PreRes)
        use  lsmrblasInterface, only : dnrm2
        IMPLICIT NONE
! INPUT
        INTEGER,INTENT(IN):: dall,maxvp
        REAL,INTENT(IN):: GVs(dall,maxvp)
        REAL,INTENT(IN):: dv(maxvp),sigmaT(dall),Tdata(dall)
! OUTPUT
        REAL:: resbst(dall)
        REAL:: res2Norm,resNSigma2Norm
        REAL:: PreRes
! PARAMETER
        REAL VsInv(maxvp),GcInv(maxvp),GsInv(maxvp)
        REAL fwdTgc(dall),fwdTgs(dall), fwdTvs(dall),fwdt(dall)
        REAL fwdTaa(dall),res(dall),resW(dall),per_res(dall)
        INTEGER jj,kk,ii
        INTEGER i,j,k
!--------------------------------------------------------------------!
! Target:
!  Calculate traveltime resiudal between observe traveltime
!  and forward traveltime which is generated by inversion result
!  using G from ref  model.
!   r=Gm-d
!--------------------------------------------------------------------!
        VsInv=dv(1:maxvp)
!--------------------------------------------------------------------!
! forward residual
        fwdTvs=MATMUL(GVs,VsInv)
        fwdt=0
        res=0
        resW=0
        per_res=0
        resbst=0
        DO i=1,dall
                resbst(i)=Tdata(i)-fwdTvs(i)
                resW(i)=resbst(i)*1/sigmaT(i) ! here, 1/sigma is for
                per_res(i)=resbst(i)/Tdata(i)
        ENDDO
        res2Norm=dnrm2(dall,resW,1)
        resNSigma2Norm=dnrm2(dall,resbst,1)
        PreRes=sum(abs(per_res(1:dall)))/dall
        write(6,*)' Mean Residual Percentage (%)',PreRes*100
        write(66,*)' Mean Residual Percentage subroutine (%)',PreRes*100

        END SUBROUTINE


        SUBROUTINE Calmodel2Norm(nar1,nar,maxvp,count3,rw,col,dv,weight,Mnorm2,MwNorm2,isTikh)
        use  lsmrblasInterface, only : dnrm2
        IMPLICIT NONE
! INPUT
        INTEGER,INTENT(in) ::nar1,nar,maxvp,count3
        REAL rw(*)
        INTEGER col(*)
        REAL,INTENT(IN):: dv(*)
        REAL,INTENT(IN):: weight
! OUTPUT
        REAL Mnorm2,MwNorm2
! PARAMETER
        INTEGER jj,kk,ii
        INTEGER i,j,k
        REAL Lm(count3),LmWeight(count3)
        INTEGER numC
        LOGICAL:: isTikh
!--------------------------------------------------------------------!
! count3: row number of L
! Formula: Lm=[a12(m1-m2)
!                              a13(m1-m3)]
! nar1: number of  nonzero value of G
! nar: number of nonzero value of Gbar-E
! dv: inversion result
! rw: value of Gbar
! iw2: row number of Gbar, not contain column
! col: column number of Gbar
! Mnorm2: 2 norm of Lm
! numC: count index of  nonzero L
! weight: in L, weight has added to Gbar, so divide weight
! Regularization Type: (T) 1st order Tikhonov ;(F) Gaussian

    IF(isTikh)THEN
        DO i=1,count3
              numC=nar1+i
              Lm(i)=rw(numC)*dv(col(numC))/weight
              LmWeight(i)=rw(numC)*dv(col(numC))
        ENDDO
        Mnorm2=dnrm2(count3,Lm,1)
        MwNorm2=dnrm2(count3,LmWeight,1)
    ELSE
          IF (((nar-nar1)/2-count3).NE. 0) THEN
               write(*,*) '  count3=',count3
               write(*,'(a,i9)') '  (nar-nar1)/1=',(nar-nar1)/2
             STOP 'L row number /= count3, in sub Calmodel2Norm!'
        ENDIF
        Lm=0
        numC=0
        DO i=1,count3
                numC=nar1+i*2-1
                Lm(i)=rw(numC)*dv(col(numC))/weight+rw(numC+1)*dv(col(numC+1))/weight
                LmWeight(i)=rw(numC)*dv(col(numC))+rw(numC+1)*dv(col(numC+1))
        ENDDO
        Mnorm2=dnrm2(count3,Lm,1)
        MwNorm2=dnrm2(count3,LmWeight,1)
        !        write(*,'(a,i19)'),'nar=',nar
        !        write(*,'(a,i19)'),'numC=',numC
        IF ((nar-numC-1).NE. 0) STOP 'nar /= numC+1, in sub Calmodel2Norm!'
    ENDIF







!        open(90,file='Lnorm.txt')
!        DO i=1,count3
!               numC=nar1+i*2-1
!            write(90,*)numC,col(numC),col(numC+1),dv(col(numC)),dv(col(numC+1)),rw(numC),rw(numC+1)
!        ENDDO
!        close(90)
!
!        open(90,file='dv_output.txt')
!        DO i=1,maxvp
!            write(90,*)dv(i),dv(i+maxvp),dv(i+maxvp*2)
!        ENDDO
!        close(90)

        END SUBROUTINE
